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In  applying  the  Wiener  multichannel  time- series  processing 
theory  to  seismic  arrays  for  teleseismic  signal  extraction,  the  transient 
nature  of  the  signal  is  not  used.  This  report  presents  the  theoretical  devel¬ 
opment  of  optimum  nontime- stationary  signal- extraction  filters  which  take 
advantage  of  the  additional  information  of  known  signal-arrival  time.  The 
matrix  inversion  is  maae  computationally  practical  by  using  an  extension 
of  the  Levinson  technique  to  determine  prediction  filters  several  points  into 
the  future.  A  sequence  of  programs  to  accomplish  the  filter  design  and 
application  is  described  anc'  illustrated  by  a  sample  problem.  Significant 
improvement  in  mean- square  error  was  not  obtained  for  this  sample  prob¬ 
lem  cor.'  jared  to  classical  Wiener  filters.  Full  evaluation  of  this  technique 
requires  additional  data  processing  and  examination  with  respect  to  criteria 
other  than  mean-square  error. 
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SECTION  I 

SUMMARY  AND  INTRODUCTION 


In  applying  the  Wiener  multichannel  time-series  processing 
theory  to  seismic  arrays  for  teleseismic  signal  extraction,  the  transient 
nature  of  the  signal  is  not  used.  This  report  presents  the  theoretical  devel¬ 
opment  of  optimum  nontime  -  statiom  ry  signal-extraction  filters  which  take 
advantage  of  the  additional  information  of  known  signal-arrival  time.  The 
matrix  inversion  is  made  computationally  practical  by  using  an  extension 
of  the  Levinson  technique  to  determine  prediction  filters  several  points  into 
the  future.  A  sequence  of  programs  to  accomplish  the  filter  design  and 
application  is  described  and  illustrated  by  a  sample  problem.  Significant 
improvement  in  mean-square  error  was  not  obtained  for  this  sample  prob¬ 
lem  compared  to  classical  Wiener  filters.  Full  evaluation  of  this  technique 
requires  additional  data  processing  and  examination  with  respect  to  criteria 
other  than  mean-square  error. 

Multichannel  data  from  an  array  of  seismometers  are  assumed 

to  be 


where 

s.  =  seismic  signal  from  seismometer  i  at  time  t 

1  j  t 

n  =  noise  from  seismometer  i  at  time  t 

lj  t 

x  =  output  of  seismometer  i  at  time  t 
19  t 

It  is  desired  to  find  the  filters  or  linear  combinations  of  t.ie  output  for  times 
t  =  1, . . . ,  k  +  m  which  are  best,  in  the  mean-square-error  sense,  for  esti¬ 
mating  the  signal  s  .  In  addition,  we  assume  that  the  signal-arrival  time 
is  known  to  be  at  time  t  =  k  +  1  for  each  seismometer. 
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The  purpose  of  this  report,  therefore,  is  to  evaluate  and 
describe  the  signal-extraction  process  developed  to  use  the  additional 
information  of  known  signal-arrival  times.  A  series  of  programs  has  been 
written  to  perform  the  processing,  and  data  recorded  at  L,ASA  have  been 
selected  and  processed.  Results  of  the  nontime -stationary  processing  have 
been  compared  with  the  usual  MCF  approach  by  theoretical  mean-square- 
error  comparisons  and  with  actual  data. 
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SECTION  II 

THEORETICAL  DEVELOPMENT 


The  notation  to  describe  the  actual  data  processing  is  more 
difficult  to  follow  than  necessary  for  an  explanation  of  the  method.  In  this 
section,  the  method  is  explained  generally;  and  the  exact  steps  of  the  com¬ 
putation  are  given  explicitly  in  Appendix  A. 

Let  the  row  vector  n^  denote  the  multichannel  pure  noise  data 

(n.  )  t  =  1,  . . .  ,  k,  which  precedes  the  signal  arrival.  The  vector  n  +  s  will 

t  2 

denote  the  multichannel  signal-plus -noise  data  (s.  +  n.  )  in  the  time  interval 

t 

t  =  k+  1,  . . . ,  k  +  m.  It  is  desired  to  find  the  least-mean- square  estimate  of  s 
by  filtering  the  data  set  (n^,  n^  +  s). 

Most  results  in  this  section  can  be  easily  seen  in  terms  of 
the  following  lemma:  Let  a  random  vector  x  with  zero-mean  be  partitioned 
in  the  form  x  =  (x  ,  x  )  with  covariance 

1  L 


n  =  Ex  x  = 


"  nn 

niz‘ 

.  n21 

fi22. 

(2-1) 


where  Ex.  x. 

i  J 

x2  are 


=  Q.j.  Given  x^,  then  the  least-mean-square  estimates  of 


T  -IT  T 

x2  =  n21  "ll  *1  *  r  X1 

(assuming  a  mear  of  0)  with  the  filters  T  defined  by 


(2-2) 
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21  11 


<r(1ll»T  -  n21T 


niir  ■  n21 


The  covariance  matrix  of  the  error  ir 


E^X2  "  X2^  ^X2  '  X2^  "  ^22  ^21  ^11  ^12  =  P 


Equations  (2-3)  and  (2-4)  are  often  expressed  in  matrix  form  as 


nn  nizi  r-rT 


n21  n22  1 


We  now  make  the  correspondence  x  =  s  and  x.  =  (n  ,  n  +  s) 
T  “  Alfa 

so  that  if  S  =  Es  s  and 


E(nj,  n2)  (nj,  nj  = 


N11  N12 


N21  N22 


with  En^  n.  =  N„,  it  follows  from  Equation  (2-3)  of  the  lemma  that  the 
optimum  filters  F  =  (F.  F.)  are  given  by 


v  •:  w- 
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VJ 


N, , 

N.  _ 

rFlTi 

o 

11 

12 

i 

N,, 

N__  +  S 

ft 

s 

L  21 

22  j 

L  2  J 

and  the  signal  estimates  are 

T  T 

s  =  F(nj,  n  +  s) 


(2-7) 


(2-8) 


with  mean-square  error  given  by 


MSE1  =  S  -  F(0  S)T  =  S  -  F2S 


(2-9) 


Although  the  problem  is  now  solved,  the  solution  is  not  in  the  best  form  for 
calculation.  We  observe  that  Equation  (2-5)  with  Xj  =  nj  and  x^  "  n 2  can  be 
modified  to 


"Nn 

Nu 

- 1 

1 

•"1 

_ i 

0 

.N21 

N22  +  S. 

i 

P  +  s 

(2  10) 


by  adding  S  to  N  and  to  P.  Then,  by  multiplying  on  the  right  by  (P  +  S)~  1  S, 
we  obtain  Equation  (2-7).  Thus,  we  have 


(P  +  s)"1  S 


(2-11) 
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so  that  the  nontime -stationary  filter  estimates  are 


/\ 

sT  =  F2(-r,  I)  (nl(  n2+  s)T 


T  T  T 

F2(n2  +sT-miT) 


(2-12) 


In  summary,  the  optimum  filters  to  apply  to  (n  ,  n  +  s)  for 

i  La 

signal  estimation  are  equivalent  to  estimating  n  from  n  and  subtracting  the 

b  1 

estimates  from  the  observed  data  n2  +  s,  then  solving  the  smaller  dimensioned 
problem  of  estimating  s  by  the  optimum  filters  operating  on  s  +  n  . 

Results  given  ir  a  later  section  in  terms  of  MSE  or  the  conven¬ 
tional  and  nontime- stationary  filters  are  more  meaningful  if  we  present  some 
theoretical  relations  which  are  known  to  hold  for  the  particular  signal  model 
used  as  the  example  in  this  report.  Comparing  MSE^  from  Equation  (2-9) 
with  conventional  Wiener  MSE  provides  an  optimistic  estimate,  for  certain 
signal  models,  of  the  expected  improvement  from  using  the  filters  F  instead 
of  the  conventional  filters.  This  comparison  is  optimistic  because  conven¬ 
tional  Wiener  MSE  values  are  derived  under  the.  assumption  that  the  data  to  be 
processed  are  (n^ +Sj,  n?  +  s2).  The  MSE2  obtained  when  the  Wiener  filters 
are  applied  to  (n^,  n2  +  s)  will  be  less  than  the  conventional  MSE.  These 
statements  can  be  made  algebraically  precise  by  again  using  the  lemma. 

Wiener  filters  f  are  conventionally  designed  as  the  optimum  filters  applied 
ton.  +  n  +  s  to  estimate  sn,  a  scalar  component  of  either  s.  or  s  . 

In  the  present  comparison,  we  restrict  s^  to  be  a  scalar  component  of  s2. 

We,  therefore,  employ  the  lemma  by  making  the  correspondence  x_  =  sn, 

C*  v/ 

X1  =  ^nl  *  si»  n2  +  8 2^  t0  °k*a*n  the  filter 

f  =  Yff1 
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from  Equation  (2-2),  where 


E  xj  xx  =  n 


<E\  V  Es0s2>  *  (V1  V  *  * 


a0  =  Ex2*2 


MSE  =  aQ2  -  fYT 


An  equivalent  problem  is  to  find  the  filter  f,  minimizing  the 


quadratic  form 


MSE  =  (1,  -f) 


d,  -f) 


Now  let 


N..  =  E  n.  n. 

i  J 


S..  =  E  s.  s. 
ij  i  J 


Then  the  Wiener  filters  f  applied  to  the  data  (n.,  r.  +  s  )  results  in  a 

^  b  U 

mean-square  error 


solsnos  strvloM  division 


mse2  =  (I,  -f)  < 


ao  Y 


T 

y  n 


n 


'21 


12 


(1.  -f) 


so  that 


MSE  a  MSE  -  (1,  -£) 
c* 


1 

0  ‘ 

S 

11 

12 

0 

21 

(1.  -O' 


If  the  signal  model  is  such  that  S12  =  0,  then  since  Sjj  is 
positive  definite,  one  has  the  heirarchy  MSE^  £  MSE2  £  MSE.  Plotted 
values  of  MSEj  and  MSE  are  discussed  in  the  next  section  for  a  specific 
example.  The  comparitive  approprite  term  MSE  was  not  evaluated. 


SECTION  III 
EVALUATION 


A.  DATA  DESCRIPTION 

Two  adjacent  noise  samples  recorded  at  LASA  subarray  B1 
were  selected  for  the  evaluation  of  the  nontime- stationary  technique. 

•  Noise  1,  25  March  1966,  04:26:12.8  to  04:34:12.7 

•  Noise  2,  25  March  1966,  04:34:13.  1  to  04:42:13.0 

Both  noise  samples  has  been  filtered  with  an  antialiasing  filter  and  resampled 
to  100  mils.  To  take  advantage  of  the  program  restrictions,  seven  channels 
of  data  recorded  at  the  center  seismometer  and  t  o  outer  rings  were  selected 
(Figure  III- 1). 

Another  antialiasing  filter  was  applied,  and  the  data  were 
resampled  to  300  mils.  Channel  1  of  noise  sample  1  was  used  to  design  a 
21-point  deconvolution  filter,  and  the  filter  was  applied  to  both  noise  samples. 
A  correlation  matrix  was  generated  using  the  300-mil  whitened  noise  sample  1. 

B.  DISCUSSION  OF  PROCESSING  RESULTS 

The  noise  correlation  matrix  was  used  as  input  data  to  the 
sequence  of  programs  described  in  Appendix  B.  Twenty-point  multichannel 
prediction  filters  were  designed  to  predic.  data  vectors  1,  2,  . . .  t  15  points 
ahead  of  the  data.  To  provide  a  better  evaluation  of  the  theoretical  quality 
of  the  multichannel  filters,  single-channel  prediction  filters  were  designed 
also.  Predictability  of  each  channel  (100-percent  minus  percent  of  mean- 
square  error)  by  the  multichannel  filters  was  plotted  along  with  the  single¬ 
channel  predictability  (Figures  III-2  and  III- 3 ).  Average  multichannel 
predictability  ranged  from  48.7  percent  at  one  point  ahead  to  16.3  percent 
at  15  points  ahead,  while  single- channel  predictability  ranged  from  12.7  to 
6.6  percent  over  the  same  number  of  points.  The  multichannel  predictability 
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decreases  rapidly  between  M  =  1  and  M  =  6;  and  from  6  points  ahead,  there 
is  a  consistent  difference  of  about  10  percent  between  multichannel  and  single¬ 
channel  predictability. 

The  single-channel  predictability  plots  in  Figure  III- 2  indicate 
that  the  21-point  deconvolution  filter  did  not  successfully  prewhiten  the  data. 
The  first  10  lags  of  the  whitened  correlation  were  plotted  and  appear  to  be 
behaving  reasonably.  It  appears  now  that  a  shorter  deconvolution  filter  would 
have  been  more  successful.  The  expense  of  returning  to  the  preprocessing 
stage  and  the  absence  of  encouraging  results  made  shorter  prediction  filters 
seem  like  the  most  efficient  way  to  overcome  this  difficulty.  Therefore, 
5-point  prediction  filters  to  predict  1  to  4  points  ahead  were  designed,  using 
the  same  prewhitened  correlations. 

Single- channel  and  multichannel  predictability  are  plotted  in 
Figure  III- 3 .  The  single-channel  predictability  on  all  channels  has  been 
significantly  reduced  by  prewhitening,  and  channel  1  is  approximately  0  as 
would  be  expected.  The  difference  between  multichannel  and  single- channel 
predictability  appears  to  remain  about  the  same. 

An  infinite-velocity  signal  model  with  gain  slop  was  designed, 
using  as  the  zero  lag  of  the  autocorrelation  three  times  the  zero  lag  of  the 
noise  autocorrelation  and  0's  for  the  other  lags.  The  signal  model  and  pre¬ 
diction  filters  were  input  to  R687,  and  10- point  multichannel  signal  extraction 
filters  were  designed  and  written  on  the  CFT  for  application.  The  values  of 
MSE  were  computed  as  described  in  Section  II  and  printed.  The  same  signal 
model  was  added  to  the  complete  noise-correlation  set,  and  a  Wiener  multi¬ 
channel  filter  to  estimate  channel  1  was  designed.  The  filter  was  designed  to 
a  length  of  35  points,  and  the  mean-square  error  was  printed. 

Since  the  data  used  in  this  experiment  has  not  been  normalized, 
the  mean-square  errors  were  divided  by  the  zero-lag  ~ignal  correlation  for 
comparison.  Plotted  were  the  normalized  MSE  from  the  35-point  Wiener 
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Figure  III- 2.  Predictability  for  20- Point  Prediction  Filters 
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Figure  III-3.  Predictability  for  5-Point  Prediction  Filters 

filter  design  (as  an  x  on  Figure  III-4)  and  the  normalized  MSE^  for  each 
point  of  channel  1  from  the  R687  output. 

Once  again,  Wiener  MCF's  were  designed  for  the  purpose 
of  comparing  mean-square  error.  Filters  of  9  points  were  computed, 
estimating  each  of  the  four  signal  points;  results  are  shown  in  Figure  IH-4. 

To  provide  a  visual  evaluation  of  the  signal-extraction  process, 
a  signal  was  added  to  each  channel  of  noise  2  at  the  same  time  interval.  The 
time-aligned  signal  plus  noise  was  input  to  R688  for  application  of  the  filters 
designed  by  R685  and  R687.  First,  10  of  the  noise -prediction  filters  were 
applied  to  the  20  points  of  noise  just  ahead  of  the  inbedded  signal  to  predict 
the  noise  in  those  10  points  of  signal  plus  noise.  The  noise  estimate  was 
subtracted.  Then  the  10-point  signal- extraction  filters  were  applied  to  the 
signal- plus -noise  prediction  error.  Figure  III- 5  shows  a  plot  of  the  seven 
channels  of  data  with  a  vertical  axis  to  indicate  the  signal  onset,  the  noise 
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SECTION  IV 
CONCLUSIONS 


There  was  no  significant  improvement  in  mean-square  error 
beyond  the  first  four  signal  points.  As  a  general  signal-processing  scheme, 
the  nontime- stationary  technique  does  not  appear  to  offer  any  significant 
advantages.  However,  the  slight  improvement  in  the  initial  points  suggests 
that  the  method  may  have  value  in  determining  first  motion. 

For  determination  of  first  motion,  the  mean-square-error 
criterion  may  have  only  limited  significance.  If  such  a  study  were  desired, 
it  would  probably  require  the  processing  of  several  events  of  lew  signal-to- 
noise  reitio. 
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APPENDIX  A 

PROGRAM  SPECIFICATION 


APPEND  JC  A 


PROGRAM  SPECIFICATION 

A.  NOISE  PREDICTION,  STEP  1 

'  The  Fr°blem  is  most  easily  solved  when  considered  as  a 
2-step  signal-extraction  process.  In  the  first  step,  the  k-sample  vectors 
of  noise  are  used  to  predict  the  noise  in  the  m- sample  vectors  of  signal 
plus  noise.  Then,  the  noise  prediction  is  subtracted  from  the  signal  plus 
noise,  leaving  signal-plus-noise  prediction  error.  In  the  second  step, 
optimum  signal-extraction  filters  are  designed  and  applied  to  the  signal- 
plus-noi-ie  prediction  error. 

To  implement  the  process,  the  first  step  is  to  design  predic¬ 
tion  filters  capable  of  predicting  noise  vectors  in  the  signal-plus-noise  area. 
To  express  this  problem  in  familiar  notation,  reorder  the  data  so  that 

-T  /-  T  — T  -  t\ 

n  =  I  n  ,  n.  1,...,n  I 
y  k-.  m  k+m-1  1  / 


N(k+m) 


the  noise  prediction  becomes  equal  to 

N12  N22  n  (2) 


If  the  set  of  prediction  filters  is  denoted  by 


then 


and 


or 


X  ■  N!2  "a 


X  N22  -  N!2 


N22X  ■  N2. 


which  can  be  written  as 


Or)  T 

Nl  '  X  =  N 


21 


From  this  matrix  equation,  the  particular  system  of  equations 
for  designing  an  NC- channel  k-point  prediction  filter  to  predict  r  points 
ahead  is  given  by 


*Nrt  N  ,  .  .  .  N  ,  " 

X, 

’n 

0  -1  -k+1 

lr 

r 

N,  N_  .  .  .  N  ,  _ 

X„ 

N  , 

1  0  -k+2 

•  • 

2,  r 

• 

= 

r+1 

• 

t  • 

•  • 

N  ,  N,  _  .  .  .  N. 

• 

• 

X, 

•  | 

• 

N 

k- 1  k-2  0  J 

L  k,  r 

.  r+k- 1 . 

A  recursive  procedure  for  the  numerical  solution  of  these 
equations  now  will  be  developed#  Using  the  Levinson  forward  and  backward 
filters  which  predict  one  point,  the  procedure  iterates  to  the  filter  predicting 
r  +  1  points  ahead  from  the  filter  predicting  r  points  ahead.  The  k-point 
Levinson  prediction  filters 
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where  P  and  P'  are  the  respective  theoretical- error  covariance  matrices, 
k  k 

Programs  to  compute  G  and  G'  from  an  efficient  procedure  already  exist. 


Using  the  notation 


r  k 

I  Xl.r 


for  the  k-point  prediction  filter  predicting  r  points  ahead,  it  follows  from  the 


k+ 1  k 

definitions  of  G  and  X  that 

r 


(A-  1) 


which  provides  a  starting  point  for  the  recursive  procedure. 
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Recursion  formulas  to  compute  Xk+J  when  Xk  is  known  will 
now  be  developed.  Assume  X  k  is  known.  Observe  that 


+  Gk  X  k 
1.  r 


holds  when  multiplying  on  the  left  by  N(k)  so  that  we  can  solve  for 


X k-  ^  _v-k  ,,  k  k 
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for  is  1,  2,  ....  k-1 


so  that 
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which  implies  that 
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for  i  =  1,  2,  . . . ,  k-1 
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Put,  using  the  above  expression  for  X  these  equations 


become 


Equations  (A-l),  (A-2),  and  (A-3)  establish  the  recursive  procedure  for 
computing  the  prediction  filters. 

The  application  of  these  filters  to  predict  the  noise  ahead  and 

the  subtraction  of  the  predicted  noise  from  the  signal  plus  noise  result  in 

signal- plus -noise  prediction  error.  The  variance- covariance  matrix  of  the 

noise-prediction  error  can  be  computed  by  observing  that  the  error  n  can 

e 

be  expressed  as 


/ 


ne  =  n(l)-n  =  n(l)-N12N22  n<2)  = 


,n(2). 


so  that  the  variance- covariance  matrix  is  given  by 
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B.  SIGNAL  EXTRACTION,  STEP  2 


The  data  from  the  array  of  seismometers  is  now  in  the  .orm 


x;  = x.  -  S = 


for  t  a  k+1,  k+2,  ....  k+m 
and  i  s  1,  2,  . . . ,  NC 

The  second  phase  of  the  process  is  to  design  and  apply  these  filters  to  the 
signal-pius-noise  prediction  error. 

The  random  vectors  are  assumed  to  be  zero  mean;  the  signal 

and  noise  are  assumed  to  be  independent.  The  variance -covariance  matrix 

of  joint  distribution  of  the  random  variables  s  and  s  +  n  is 
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where  the  only  difference  between  this  equation  and  the  usual  Wiener  equation 
is  that  is  not  Toeplitz  in  form  and  a  general  matrix  inversion  will  replace 
the  Levinson  iteration. 

Two  approaches  to  multichannel  least-squares  filtering  of 
seismic  data  have  been  presented  here.  The  signal  estimate  can  be  com¬ 
puted  by  using  the  general  solution  F  as  follows: 

k+m 

r  =  V  f-T  X  ... 

P  ^  jp  k+m-j+1 
j=l 

or  equivalently  from  the  noise  prediction  solution  by 
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and 
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APPENDIX  B 

PROGRAM  DESCRIPTIONS 


The  nontime -stationary  signal- extraction  procedure  is 
implemented  by  a  series  of  three  computer  programs  developed  for  use 
on  the  IBM  7044  processor.  The  first  program,  R685,  designs  the  pre- 
dection  filters;  the  second  program,  R687,  designs  signal- extraction  filters 
using  signal-plus-noise  prediction  error;  the  last  program,  R688,  applies 
the  filters  developed  in  the  first  two  programs  to  data.  A  more  detailed 
description  of  the  capabilities  and  limitations  of  these  programs  follows. 

A.  R685  -  PREDICTION  FILTER  DESIGN 

A  noise  correlation  set  is  input,  and  prediction  filters  are 
designed.  For  . 

NC  =  number  of  channels 
NFP  =  number  of  filter  points 
M  =  number  of  points  to  be  predicted  ahead 
L  -  number  of  Jags  in  correlation 

the  following  restriction  exists: 

NC2  *  M  5  2000 


Since 


L  a  M  +  NFP 


the  size  of  the  filter  and  the  number  of  filters  designed  are  restricted.  Of 
course,  the  restriction  can  be  relaxed  by  a  program  modification  if  longer 
filters  or  more  distant  predictions  appear  to  have  merit. 


The  program  provides  a  visual  evaluation  of  the  predictability 
of  the  data  by  generating,  for  each  channel,  a  plot  of  the  predictability  per¬ 
centage  (Figure  III-3). 

B.  R687  -  SIGNAL  EXTRACTION  FILTER  DESIGN 

The  prediction  filters  designed  by  R685  and  a  signal  model 
are  input,  and  signal  extraction  filters  are  computed.  For 

M  =  number  of  filter  points 
NC  =  number  of  channels 

the  following  restriction  exists: 

(NC  *  M)2  <;  5500 

The  program  also  computes  the  diagonal  elements  of  the  variance-covariance 
matrix.  Therefore, 

MSE,  =  S  -  S(S  +  N  )_1  S 
1  e 

for  evaluation  of  the  filters'  theoretical  effectiveness. 

C.  R688  -  FILTER  APPLICATION 

The  prediction  filters,  signal- extraction  filters,  and  data 
se;  are  input;  and  the  filters  are  applied.  The  prediction  filters  are  applied 
to  the  first  data  vectors,  and  the  predicted  noise  is  subtracted  from  each 
of  the  next  vectors  of  signal  plufi  noise.  Then,  the  signal- extraction  filters 
aro  applied  to  the  signal-plus -noise  prediction  error.  The  input  data, 
noise  estimate,  signal-plus-noise  prediction  error,  and  extracted  signal 
then  are  plotted  on  the  same  set  of  axes. 
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